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Abstract. We report on a new study of the velocity distribution in N-body simulations. We investigate the center-of-mass and 
internal kinetic energies of coarsening cells as a function of time, cell size and cell mass. By using self-similar cosmological 
models, we are able to derive theoretical predictions for comparison and to assess the influence of finite-size and resolution 
effects. The most interesting result is the discovery of a polytropic-like relationship between the average velocity dispersion 
(internal kinetic energy) and the mass density in an intermediate range of densities, /ej nt oc g 2 ^. The exponent n measures the 
deviations from the virial prediction, ^ v iriai = 0. For self-similar models, n depends only on the spectral index of the initial 
power spectrum. We also study CDM models and confirm a previous result that the same polytropic-like dependence exists, 
with a time and coarsening length dependent 77. The dependence K mi (g) is an important input for a recently proposed theoretical 
model of cosmological structure formation which improves over the standard dust model (pressureless fluid) by regularizing 
the density singularities. 
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1. Introduction 

The use of N-body simulations has proved a useful tool in 
the investigation of both the cosmological structure formation 
and the evolution by self-gravity. The main interest has been 
concentrated in properties of the spatial distribution of mat- 
ter (mass correlations, void distribution, morphological fea- 
tures...), while the kinetic properties have received compara- 
tively little attention, without doubt due to the larger difficul- 
ties to obtain reliable ki netic meas urements from real data with 
which to compare ( Watk ins et all2 002). 

The kinetic measurements addressed till now with N- 
body simulations have pertained the quasili near velocity field 
(see e.g. the reviews by iDekell Il994l and iBernardeau et alJ 
2002t an application c losely related to the present work is 
Seto & Sugivama 2001), the pairwise relati ve velocity (jP eebles! 
1980; fo r recent applications, see e.g. IStrauss et alJ [1998; 
Feld man et alJ l2003h. and t he velocity dispersion of halos 
(e.g. iKnebe & Mulled Il999l in connection with the present 
work). Our work addresses the center-of-mass velocity of 
coarsening cells (macroscopic kinetic energy), as well as the 
velocity dispersion of the particles inside the cell (internal 
kinetic energy). The coarsening cells are randomly centered 
and of variable size (probing both the linear and the non- 
linear regimes); in this way our analysis does not suffer the 
arbitrariness intrinsic to the definition of clusters and halos 



jKlvpin & Holtzmannlll99'fl and refs. therein), and essentially 
all the simulation particles are employed in the determina- 
tion of the quantities. We are aware of two works where a 
simil ar analysis of N-body simulations ha s been performed 
jKepner et alJll997t iNagamine et alJbOOlh . motivated differ- 
ently than ours. The connection to our work is explained in 
Sec.|5]in detail. 

The present study focuses on the dependence of the cell ki- 
netic energies on the cell mass density. The main motivation is 
the application to models of cosmological structure formation 
by self-gravity. The most widel y used theoretical model is the 
dust m odel (pressureless fluid) JPeebleslll980l Ipadmanabhan 
11994 . whic h has been studied intensively (see e.g . the re- 
views ISahni & Colesl fl995l IBernardeau et all I2002) but has 
the sho rtcoming of producing singularities. Some recent 
works teuchert & Domfnguezlll998|: lAdler & Bucher3 | 19991: 
" iMaartens et allll999t IPomuiguezlkOOO ; 



Bucherte talJll99 . _ ..... _ 

Morita & Tatekawa l200lt lDomfnguezll2002t iTatekawa et all 
2002) have proposed a novel approach. O ne of its features is the 
abilit y to derive adhesion-l i ke models jKofman & Shandarir 
198& iGurbatov et alJll989t iKofman et alJll992t iMelott et al 
1994; Sathvaprakash et al. 1995), and to offer a possible expla- 
nation of the physical origin for the "adhesive" behavior which 
regularizes the mass density singularities of the dust model. In 
these improved models, the internal kinetic energy brings about 
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the "adhesive" behavior provided it can be approximated as a 
function of density and/or the gradients of the velocity field. 

Another goal of thi s work is to confirm the results by 
Dommguezl ill 9991 fcoCBi. where a polytropic-like dependence 
between internal kinetic energy and mass density is found in 
AP 3 M simulations of CDM models. We indeed corroborate this 
finding in PM simulations of CDM models and also of self- 
similar models. The latter are particularly amenable to a theo- 
retical analysis and allow the identification of the influence of 
finite-size and resolution effects in the measurements. We con- 
clude that the polytropic-like dependence is unlikely to be an 
artifact of the simulations. 

The paper is organized as follows: in Sec.[2]we work out 
the theoretical predictions for the density dependence of the 
macroscopic and internal kinetic energies. In Sec. [3] we de- 
scribe the simulations and the method how we measure the ki- 
netic energies. In Sec.0]we present the results of the analysis. 
Sec.|5]contains a discussion of the results and the conclusions. 

2. Theoretical background 

Let a(t) denote the cosmological expansion factor, m the mass 
of a particle, and x a and u a the comoving position and peculiar 
velocity, respectively, of the a-th particle. W{-) is a (normal- 
ized) smoothing window. Then, given a comoving smoothing 
scale L, the coarse-grained mass density field, velocity field 
and density of internal peculiar kinetic energy are defined re- 
spectively as follows: 



g(x, t; L) = 



gu(x, t; L) 



[a(t)L\ 



[a(t)L] 



a=l 



■ J] Ua(t) W 



X a (t) 



Qr=l 



' x„(?) 



(1) 



K int (x, f, L) = 



2[a(t)L] : 



J]\u a (t)-u(x,t;L)\ 2 W 



■ x ff (f) 



The density of macroscopic 1 peculiar kinetic energy is defined 
as ^mac = (l/2)p|u| 2 , so that one can write K l0 , = + K mac , 
where K wt is the density of total peculiar kinetic energy: this 
is given by the same expression as K mt above but dropping u. 
Whether K mlLC or K mt dominates the contribution to K tot means 
respectively that the particle velocities u a are controlled by 
the center-of-mass motion of the coarsening cell as a whole 
or by "internal" motions within the cell (th e ratio K mac IK\ nt 
was in troduced as the cosmic Mach number by Ostriker & Suto 
Jl990l) V Indeed, K mt is a much better probe of the dynamics at 
the small scales than K, mc , because in the latter there can be 
extensive cancellations in the vectorial sum defining u. 

The purpose of the present study is the relationship between 
Ki nt and the field g. Our main interest is the density dependence 
of (Ki nl ) g , that is, the average of K mi conditioned to a given 
value g of the coarse-grained mass density, 

(K int (x,t;L)6 (?) (g(x,t;L) - g)) 



K int (g,t;L) := (K int ) e = 



(6^(g(x,t;L)-g)) 



(2) 



where (•} denotes ensemble average, which by translational in- 
variance must be x-independent, and £ (3) (-) is Dirac's delta 
function. In the same way is defined the conditioned average 
Kmac := (K m!1 c)e- Since Ki nt receives contributions from the 
highly non-linear regime, the density dependence cannot be 
computed analytically in general. Nevertheless, a lot can be 
learned by way of suitable approximations, whose validity will 
be checked by comparing with simulations. 

To simplify the theoretical discussion, we consider a self- 
similar cosmological model: an Einstein-de Sitter background 
and an initial Gaussian distributed density field with power 
spectrum P(k, t m if m i) = Ak", with the bounds n > -3 (so 
that density fluctuations do not receive a divergent contribution 
from k — > 0) and n < 4 (imp osed by the unavoidable grainines s 
due to the point particles) (Peebles 1980; Padmanabha nl 1 995t) . 
The conclusions should apply qualitatively unaltered to a more 
realistic case. Let cr 2 (L, t) - (5 2 (x, t; L)) denote the variance of 
the density contrast smoothed on the scale L (6 :- (g/gb) - 1 is 
the density contrast, with gb :- (g) the background density), 



cr 2 (L, t) 



r dv 

J (2^r 



dk 



P(k,t)\W(Lk)\ 2 . 



(A tilde will denote Fourier transform). A (comoving) scale of 
nonlinearity, r n \(f), is defined by the condition 



r 



dkP(kj) = 1. 



k<2ji/r„i 

In the self-similar model, the physical properties which do 
not depend on the short and larg e distance cutoffs of P(k) ex- 
hibit a simple scaling behavior (Peebles 1980; Padmanabhan 
11994 . The only physically relevant parameters on which they 
can depend are gb(to) (through the cosmological background), 
a'ni(fo) (through the initial conditions), and the gravitational 
constant G (through the dynamics); self-similarity means that 
the time to is arbitrary. Thus, in combination with dimensional 
analysis 2 , it is shown that cr(L, t) is a function of the sin- 
gle quantity L/r n \(t); in the linear regime, L » r n \, one has 
cr 2 oc (L/r n i) - " -3 . Similarly, K[ nt / mac can be written in the fol- 
lowing suggestive form: 



(3) 



Kint/macfe t\ L) = gb(dr n \) Fint/mac | 8, 

r nl 



Here F( -) is a dimensionless function of its dimensionless argu- 
ments, and gb(ar n \) 2 is of the order of the Hubble-flow kinetic 
energy in balls of (comoving) radius r a \. The unconditioned av- 
erages (/fint) = (-K"int) and (A" mac ) follow the same scaling but 
without the 5-dependence. Deviations from this scaling behav- 
ior would mean a dependence on extra variables, e.g., on short 
or large distance cutoffs in P(k) or in the dynamics (as occurs 
with, e.g., numerical simulations). 

The task now is to characterize the functions Fj nt / mac (-). Let 
A and R denote the short and the large distance cutoffs, respec- 
tively, so that we take A <K r a \,L <sc R. Roughly speaking, we 
can say that K mac is determined by the motion at scales between 



Borrowing the nomenclature from statistical physics. 



2 Here it is to be distinguished between the dimensions of comoving 
and physical lengths, differing by a factor a(t). 
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infrared (IR) cutoff R. For the macroscopic kinetic energy we 
then estimate 



n >-l 



* ma c ~ g{\u\\ ~ g b {aRY(T l (R)(\ + 6), 
R <#mac> ~ Q b (aR) 2 cr 2 (R). 



(4) 
(5) 





r nl 

1 














K int 




K mac 



Fig. 1. Sketch representing the relative contribution of the dif- 
ferent length scales to /c mac and Kj nt . The scales A, R and r n \ are 
intrinsic to the system, the movable length L is the observation 
resolution. 

L and R, whereas K m is dominated by the motion at scales be- 
tween A and L. 

The hierarchical, bottom-up scenario exhibits a monoton- 
ically growing length scale, r n \{t), which is roughly propor- 
tional to the size of the largest collapsed clusters at time t. The 
bottom-up growth of structure by self-gravity can be sketched 
in the following picture: particles get trapped in clusters so that 
(i) the evolution above the cluster scale is dominantly ruled by 
the motion of each cluster as a whole ( = "effective particles") 
in the gravitational field of the other clusters, and (ii) the evo- 
lution below the cluster scale is driven mainly by the scales 
~ cluster size. This means a dynamical decoupling between 
scales above and belo w the cluster size, ~ r n \ (this idea has 
been explored by Domfnguez ( 2000, 2002) in order to improve 
the models of structure formation), and implies that the short- 
distance cutoff A is irrelevant. Depending on which of the two 
motions, (i) or (ii), contributes mostly to the particle velocity, 
there arise several possibilities (see Fig.Q: 

— If —3 < « < —1, it will be argued that the largest scales 
dominate. Thus, /c mac should be determined by the scales 
~ R and K\ nl by the scales ~ L. 

- If -1 < n < 4, the scales around r n \ provide the prevalent 
contribution because A is irrelevant and there is no other 
privileged scale between A and r a \. We can distinguish in 
turn two cases: 

- Linear regime, r n \ <K L: the value of /e mac is mainly set 
by scales ~ L, and that of A-; nt by scales ~ r n \. 

- Nonlinear regime, L <K r n \: now /c mac is dominated by 
scales ~ r„i, and K mt by scales ~ L. 

This description will be now elaborated in somewhat more 
detail. 

Case I: — 3 < n < — 1. There is so much power initially at the 
large scales, that the contribution of the linear modes (r n \k <sz 1) 
to the variance of the macroscopic velocity diverges if n < - 1 : 



<|u| 2 ><* 



c/k k" 2 + nonlinear contribution. 



Hence, both u and the particle velocities u ff will be mainly de- 
termined by the modes around the largest available scale, the 



For the internal kinetic energy, it proves useful to consider sep- 
arately the linear and nonlinear regimes: 

- Linear regime, L » r a \. Our hypothesis is that the main 
contribution to K mt comes from scales ~ L, so that we can 
employ the linear solutions to the gravitational instability 
to evaluate Eqs. 0. The calculations are collected in the 
Appendix; the final result reads 



*int - (^int) 



1+B1- 



6 l 



crHL) 



(6) 



In this expression, (K mt ) is the unconditioned average in- 
ternal peculiar kinetic energy in the linear approximation, 
Eq. iA.2\ . B is a dimensionless negative numerical coef- 
ficient, Eq. ( IA.3I . Both (Ki nt ) and B are IR-convergent be- 
cause the /^-dependence cancels from the difference u a -u, 
so that Eq. (0 holds. Indeed, we derive the scaling behavior 



<^int> « QbiaLfaHL). 



(7) 



Nonlinear regime, L <K r n \. In high density coarsening 
cells, 6 » 1, the Hubble flow is negligible and pecu- 
liar velocities are approximately equal to the physical ve- 
locities. Assuming stationarity, the conditioned differences 
(|Uo,-u| 2 ) e can be expected to be given by local virialization 
on scales ~ L, so that 



3|2 



1 G[g(aLy] 



{aLy 



aL 



QbiaLfd 2 , 



(<5» 1). 



(8) 



The opposite limit of low density, < 1 + 6 <K 1, requires 
a model for the expansion of voids. The simplest model 
would set |u ff - u| ~ hL typically, with h » a quantifying 
the void expansion speed. Then 



Q b (aL) 2 (\ +5), (0< 1 +(5<s 1), 



(9) 



predicting a low velocity dispersion. 
To compute the unconditioned average (K mt ), it must be 
noticed that, although the internal kinetic energy of high- 
density cells is very large, the number of low-density cells 
is much larger, and it is not clear which of the two compet- 
ing effects dominates. If one would assume that the main 
contribution comes from high-density cells, then Eq. 
would yield 



(K int ) ~ Q b (aL) 2 cr 2 (L). 



(10) 



Case II: -1 < n < 4. The linear modes do not lead to diver- 
gences and the velocities are now determined mainly by the 
scale r n \. 
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- Linear regime, L » r m . Now the linear solution states that 
the smoothed velocity u will be determined by scales ~ L, 
i.e., the smallest scales it probes, and thus one obtains 

(ID 

On the other hand, our hypothesis is that the particle veloc- 
ities u a are controlled by scales ~ r n \. Hence, we write the 
estimate 

2 Gg b (ar m f 2 

(\u a \ } e {ar a \Y, 

ar m 

assuming that the scales ~ r B \ are rather insensitive to the 
small density fluctuations at the scale L » r m . Then 

Kim ~ Qb{ar n \f(\ + $), (K int ) ~ g h (ar nl f. (12) 

- Nonlinear regime, L <K r n \. The conditioned difference 
(|u ff -u| 2 ) e is controlled by the local dynamics on scales ~ L 
like in the previous case, and expressions dHl-llOi are valid 
now too. The motion of the coarsening cells as a whole, 
and thus (|u| 2 ) e , is however assumed to be determined by 
the large scales ~ r n \, so that 

Kmac ~ Qb{ar D \f(l + 6), (K mac ) ~ Qbiarri) 2 , (13) 
in analogy to Eqs. (I12> . 

3. Simulation and analysis method 

The simulations of the self-simila r models a re described in full 
detail elsewhere (Melott & Shan darinlll993l) . They consist of 
a cubic box (periodic boundary conditions) of comoving side- 
length R containing N = 128 3 particles. The dynamical evo- 
lution was computed using a PM algorithm on a grid with 
Nyquist wavenumber = URn/R. The background cosmolog- 
ical expansion followed the Einstein-de Sitter solution and the 
initial conditions were generated by using the Zel'dovich ap- 
proximation for a Gaussian random field with a scale-invariant 
power spectrum P(k) oc In the present work, the values 
n = -2,0, +1 were considered and the data were studied at 
three different times, corresponding to a scale of nonlinearity 
r m = R/64, r n \ = R/16, and r m = R/4, respectively. In each 
case four independent realizations of the initial conditions were 
evolved. 

Two CDM models were also addressed with a PM algo- 
rithm: flat CDM (Q m = 1.0, Q A = 0.0, o- 8 * 1), and open 
CDM (Q m = 0.3, Q A = 0.0, cr 8 * 1). Each simulation con- 
tained N = 256 3 , and for each model two different box-sizes 
were considered, R - 128 Mpc and R — 512 Mpc. 

Starting from the coordinates {x a , u a } a= i...N provided by the 
simulation, the definitions Q were applied with a cubic top- 
hat window, 

W(z) = 0(1 - 2|zi|) 0(1 - 2\ Z2 \) 6(1 - 2|z 3 |), (14) 

where #(•) is the step function. The coarsening was imple- 
mented efficiently by covering the simulation box with a cubic 
grid; in order to probe a minimum of 2048 coarsening cells, the 
grid was shifted randomly, when needed. It was checked that 



the results are insensitive to the use of spherical top-hat win- 
dows instead. A comparison with the use of a Gaussian win- 
dow was also carried out. We observed a difference only for 
the case n = +1 at the earliest time, r m = R/64, because then 
a 2 diverges for a top-hat window (in the simulation, the ar- 
tificial point-particle discreteness regularizes the singularity). 
Hence, for this single case we disregarded the use of the cu- 
bic top-hat window as unphysical and used instead a Gaussian 
window (which is computationally much more costly). 

The explored values of the coarsening length L were 
equally separated in a logarithmic scale and they ranged from a 
maximum (R/5) down to a minimum « 1/2, where t :- R/N 1 ^ 3 
is the average interparticle distance. This results in scatter plots 
"kinetic energy vs. density". To compute the constrained av- 
erage (|2j, the data for the kinetic energy were binned into 
40 subintervals according to the value of the density g; bins 
containing less than 10 data points were disregarded. It was 
checked that the conclusions are robust against the amount of 
binning by varying the number of bins. The constrained aver- 
age was identified with the mean of each bin. The amount of 
scatter about this (global) mean is represented in the plots by 
scatter bars, which extend between the mean of those kinetic 
energies above the global mean, and the mean of those kinetic 
energies below the global mean. We find that this method rep- 
resents the scatter of the data in the log-plots more faithfully 
than the estimation through the variance. 

We checked the algorithm in various ways. It was applied 
to an ideal gas simulation: the results for the dependence of 
Ki nt on the density agreed with the ideal gas equation of state. 
Another check was to restrict the coarsening procedure to a 
subvolume of the simulation box (1/64 of the total volume) 
for some sample cases. As expected, we find that the data are 
somewhat noisier because of the reduced number of particles, 
but the conclusions remain the same. 

4. Results 

There are a couple of general remarks which hold for all the 
Subsecs. to follow. First, the small length scale ( enters in the 
results via mass-resolution and force-resolution effects: the 
first effect refers to the presence of a minimum non-vanishing 
mass — that of a single particle. This affects the computation 
of the unconstrained averages (K mt / m . dc ) due to undersampling 
of the cells with a mass smaller than this minimum, which also 
sets a lower bound on the value of the density at a given fixed 
L when computing the constrained averages, K[ nt , K mac . In par- 
ticular, it renders all results concerning the nonlinear regime 
(L < r„i) at the earliest probed time (r n i = R/64) rather unre- 
liable, since then r n \ « t. The second effect, force resolution, 
implies that the relative force over two particles decreases when 
they are closer than the mesh spacing of the PM algorithm, * (, 
and the velocity dispersion below this scale does not grow as 
much as it would if ( — > 0. All in all, these two effects tend to 
artificially reduce the value of the kinetic energy, in particular 
of /Tint, being more sensitive to the small scales. The theoret- 
ical discussion in Sec. [2] suggests this effect to be particularly 
noticeable when n > -1 and at the earliest times, as will be 
indeed observed. 
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Second, the influence of the "cosmic variance", i.e., of the 
fluctuations in measured quantities from one realization to an- 
other, is the strongest when n = -2. For clarity, however, we 
will show in the plots the results of a single realization, since 
the other ones yield almost identical results. 

4.1. Self-similar models: mass density variance 

For reference purposes, Fig. [5] shows the measurements of 
cr 2 (L). The results collapse well on a single function of L/r n \. 
At the earliest time and the smallest lengths, one can ob- 
serve the beginning of the crossover to the Poissonian behavior, 
cr 2 {L) cc ZT 3 , induced by the small-scale discreteness. At large 
L, one recovers the linear scaling, cr 2 (L) cc Zr"~ 3 ; due to finite- 
size effects, the case n — -2 exhibits a slight departure away 
from this dependence. 

4.2. Self-similar models: the unconstrained averages 

The average (^ ma c) is the kinetic equivalent of cr 2 . In Fig. [3 
we observe that the data for (K m - dc ) do not follow at all the scal- 
ing behavior (0 when n - -2. As explained in Sec. [2] this is 
due to finite-size effects: we have checked that the data follow 
instead the dependence (|5}- The data for the other two cases, 
n — 0, +1 on the contrary, obey the expected scaling, Eq. dl It 
when L » r n \, and Eq. (1131 when L <k r n \. A major departure 
in the three cases is observed at the earliest time (r„i = R/64-) in 
the nonlinear regime (L < r n \) due to the undersampling prob- 
lem mentioned above. This was confirmed by artificially re- 
moving from the estimate of the averages those cells with less 
than a given number of particles, yielding the same behavior in 
(^mac) as detected in the plots. 

As remarked in Sec. |2] the average (Ki Bt ) is more sensitive 
to the small-scale dynamics than (K mac ) or cr 2 are. This average 
also suffers the same undersampling problem as (K m!lc ). But 
the resolution effects are somewhat larger and prevent the data 
from following the scaling (0 perfectly. One can recognize a 
tendency for these effects to become less relevant in time, and 
to be more important for larger values of the spectral index n, 
in agreement with the theoretical discussion. 

4.3. Self-similar models: the constrained averages 
4.3.1 . Check for self-similarity 

We first checked if the measured K mt / lmc followed the self- 
similar scaling relationship Q- The conclusions are almost 
the same as derived above with the unconstrained averages 
(^int/mac): Kmac follows self-similarity very well (except if n — 
-2), while A:i n t follows it a bit less well. The important dif- 
ference is that departures from self-similarity are (sometimes 
substantially) smaller than in (K ml / m . dc ), see Fig.|4] The reason 
is that Kint/mac suffer the undersampling problem due to finite 
mass-resolution only in the small-<S end of each curve or at the 
earliest time. When the number of particles in the cell is large 
enough, force-resolution is likely the main effect and it does 
not appear to spoil self-similarity so much. 




Fig. 2. Abscissa: log(L/r m ). Ordinate: log<x 2 (L) at the three 
probed times. Solid line: expected functional dependence in the 
linear regime, cr 2 oc Zr"~ 3 . 



4.3.2. K mac 

Fig. [5] shows the function K m . dc (S) at the different times and 
coarsening lengths probed for the spectral index n — -2. The 
scaling Eq. (0} is obeyed well at all times, although we find 
small fluctuations around the linear dependence between real- 
izations. This is due to the strong dependence on the IR cutoff. 
It is also responsible for the lack of collapse of the three plots 
on a single function, the deviations following indeed the law 
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Fig. 4. Abscissa: log(l + 6). Ordinate: log[K mt /gb(ar n i) 2 ]. The 
plots are arbitrarily shifted in horizontal direction for clar- 
ity. From left to right, L/r nl * 2.67, 0.40, 0.08. When self- 
similarity holds, Eq. (0, curves corresponding to the same ratio 
L/ r n i fall on top of each other. 

regime. The theoretical scaling Eq. Jl 3i is also very well fol- 
lowed within scatter bars. For the largest values of L considered 
in the plot, a slight tendency is noticeable away from the theo- 
retical dependence, which is more obvious for n — + 1 . 



-4 1 1 1 * 1 1 

-2 -1.5 -1 -0.5 0.5 1 

Fig. 3. Abscissa: log(L/r„i). Ordinate: \og[{K) I Qb(ar n \) 2 ] at the 
three probed times: the lower points in each plot correspond to 
{K mi ); the upper points represent 10 3 (K m . dc ). The solid line is 
proportional to L 2 cr 2 (L), Eqs. J7I10I1 li . 



R 2 cr 2 (R) of the factor in Eq. 10}. The slight deviations around 
(5 = 0, most noticeable at the earliest time, correspond to the 
largest values of L and are likely due to finite-size effects too. 

Fig-Elshows K mdc {6) in the linear regime when n — 0, +1. 
The behavior ( II Q is obeyed very well within the scatter bars; 
in the case n = 0, a systematic trend away from the expected 
data collapse is observed. Fig. [71 corresponds to the nonlinear 



4.3.3. A-jnt, linear regime 

Fig.[8]shows K mi {6) for n — -2 in the linear regime. Deviations 
from the theoretical prediction, Eq. with {K mt ) and B com- 
puted with Eqs. JA.2I IA.31 . are noticeable. The probable rea- 
son for the discrepancy is that the linear regime is not really 
being probed: one can observe a large asymmetry 6 — > -6 in 
the plot, in contradiction with a centered Gaussian distribution 
for 6. (That deviations from Gaussianity are indeed important 
has been checked by estimating the probability distribution of 
6 from the simulations). The most linear case we probe for 
n = —2 (corresponding to r n \ = R/64 and L = R/5) yields 
cr ~ 0.2, which satisfies the asymptotic condition cr <K 1 only 
marginally. Thus, the ultimate origin of the problem is the IR 
cutoff imposed by the simulation box which limits the maxi- 
mum value of L. 
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Fig. 5. Abscissa: log(l+(5). Ordinate: log[K m!ic /gb(ar n i) 2 ]. Solid 
line: theoretical functional dependence, Eq. ©. All the probed 
coarsening lengths, at the three probed times (increasing from 
top to bottom) are plotted. The scatter bars are barely depen- 
dent on time, coarsening length and density; for clarity, they 
are shown only at the earliest time. 

Fig- represents K m i(5) when n = 0, +1. According to 
Eq. J12l> . Ki nt is determined by the velocities at scales r n \, which 
is ~ i at the earliest time (corresponding to the plotted data). 
Thus, the results are not too reliable in principle, as discussed 
above. Indeed, the collapse of the data on a master curve is 
marginal within the scatter bars, K mt exhibits a trend to decrease 
with decreasing coarsening length — consistent with the arti- 
ficial reduction of kinetic energy at small scales by resolution 
effects. The linear dependence predicted by Eq. fl!2i is not ob- 
served, a curvature being evident. But we cannot conclude if 
this is due to resolution effects or because the derivation of 
Eq. d!2i relies on too simple arguments. 



1.4e-02 



1.2e-02 - 



1.0e-02 



i.0e-03 



6.0e-03 - 



7e-03 



5e-03 



3e-03 



0.97 




1.03 



Fig. 6. Abscissa: 1 +6. Ordinate: K mac /gb[aLcr(L)] 2 . Only plot- 
ted are the coarsening lengths such that cr(L) < 0.05 if n = 
and cr(L) < 0.03 if n = +1. The theoretical functional depen- 
dence, Eq. i fTTV is a constant. The scatter bars extend beyond 
the plotted area. 



4.3.4. A- int , nonlinear regime 

Fig. shows K mt in the nonlinear regime at the latest time, 
Ail = R/4. For each coarsening length, the theoretical func- 
tional dependence, Eq. (jHJl, is obeyed only at the largest den- 
sities, if at all. At intermediate densities 3 , Fig. II ll shows that 
the data can be better fitted by the following polytropic-like 
behavior: 

K inl ~ eh (aL) 2 a 2 (L)6 2 -\ (15) 

where the exponent rj is detectably different from the virial pre- 
diction {rj = 0). A fit by eye yields the rough values r\(n = 
-2) * 0.3, i](n = 0) ~ 0.5, rj(n = +1) * 0.6 — the amount of 
scatter does not allow a more precise determination of rj; if the 
scatter is forgotten, it could even be that rj depends slightly on 
L, as seems to be the case when n = -2, Fig. 1101 

At the intermediate time, r n \ = R/16, a similar behavior 
can be detected, although the polytropic-like dependence d!5l > 
can be discerned only with some difficulty, Fig.^1 Apparently, 

3 The data tend to follow the fit also for very small densities, but 
mass-resolution effects are large then. 



the dynamical evolution has not proceeded so far that the inter- 
mediate regime with rj + can be clearly detected without 
interference of finite-mass effects: they show up in this case as 
an artificial reduction of K ml when 6 is small enough, i.e., when 
there are only a few simulation particles in the coarsening cells, 
see Fig.fT2l 

4.4. CDM models 

Regarding the CDM models, we have investigated only the re- 
lationship K mt (5) in the nonlinear regime at the present epoch, in 
order to assess the possibility of a polytropic-like dependence 
as in the self-similar models. Fig^]shows the measured Ki at as 
a function of the density for the flat CDM model simulated in 
a box of sidelength 128 Mpc. Now there is no reason to expect 
an exact scaling behavior like (0 and the plots cannot be made 
to collapse on a single function. Nevertheless, a dependence 
fint K 6 2 ri fits well the data for 5 large enough, with a scale- 
dependent exponent, rj(L), which decreases with decreasing L 
and ranges between and « 0.5 for the lengths L which we 
probed (L ^ 0.5 Mpc, see Sec.[3jl. 
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Fig. 7. Abscissa: log(l +6). Ordinate: log[/c mac /£>fc(arni) 2 ]- Omac 
for the case n — is shifted by a factor 10 3 ). Solid line: the- 
oretical functional dependence, Eq. Jl 31 . Only plotted are the 
coarsening lengths such that cr(L) > VlO. The scatter bars are 
barely dependent on time, coarsening length and density; for 
clarity, they are shown only for n — 0. 




Fig. 8. Abscissa: 6/cr(L). Ordinate: Ki nt /gt,[aLcr(L)] 2 . Solid 
line: theoretical prediction, Eq. (J6J. Only coarsening lengths 
are plotted such that cr(L) < 0.3. 



Fig. 1141 shows the measurements of K; nt of the flat CDM 
model in two different simulation boxes (128 and 512 Mpc 
sidelength, respectively). The worse mass-resolution in the 
largest box implies that, for a given L in the highly non-linear 
regime, the minimum measurable value of 6 is larger. It also 
means that the absolute value of K mt is smaller. The interesting 
finding is that, if K ml is multiplied by a factor * 10 3 , the plots 
corresponding to different simulation boxes but to the same 
coarsening length superpose each other. 

The conclusions extracted from the open CDM models 
are qualitatively identical to those reached with the flat CDM 
model, and the numerical values for the exponent 77 are very 
similar. 



2.4e-02 



2.0e-02 



1.6e-02 



1.2e-02 



n=+0 





7e-03 



1.03 



Fig. 9. Abscissa: \ +6. Ordinate: K mt /gb(ar n {) 2 . Theoretical pre- 
diction: a linear dependence. Only coarsening lengths such that 
cr(L) < 0.05 if n = and cr(L) < 0.03 if n = +1. 

5. Discussion and conclusions 

In the previous Sec. we have measured the macroscopic and the 
internal kinetic energies of cubic cells as a function of time, cell 
size and cell mass for different cosmological models. The use 
of self-similar models simplifies the task of comparing with 
theoretical results, when available. In particular, the scaling re- 
lationship (|3j is useful to assess unphysical dependences on the 
unavoidable additional length scales introduced by the simula- 
tion procedure, namely the box sidelength, R, and the mean 
interparticle separation, £, It must be noticed that, even when 
the results obey the scaling (0, this does not imply irrelevance 
of these extra length scales: one can conclude at most that R 
and i could enter in the result solely as the combination R/{, 
or equivalently, as N, the total particle number. We have not 
explored explicitly the influence of such a dependence on N. 
Nevertheless, from our results we can obtain some hints about 
how well they reproduce the limit N — > +00. 

We first considered the unconstrained averages, (K mt / mac ). 
We find that (K mdc ) for n — -2 suffers from strong finite-size 
effects in a predictable manner. The mean internal kinetic en- 
ergy, {K mi ), however, is strongly affected by resolution effects, 
the more so the larger « is, and its measurement is therefore 
unreliable. Next we considered the constrained averages, n ml , 
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Fig. 10. Abscissa: log(l +6). Ordinate: log[K im /g b (aL) 2 ]. Solid 
line: theoretical functional dependence, Eq. (|8J. Plotted only 
coarsening lengths at the latest time, r a \ = R/4, such that 
cr(L) > 3. For clarity, the scatter bars are omitted; see Fig.lTTl 



Kmac- In general, they are much less affected by resolution ef- 
fects, which are "localized" to very small mass densities or to 
the earliest time, being more conspicuous for n — + 1 . 

The macroscopic kinetic energy, «inac> of the case n — -2 
depends strongly on R, as predicted theoretically. However, this 
does not break self-similarit y of the amplitude of densit y fluc- 
tuations, as shown bv lTain & Bertschin ger ( 19961 Il998l) . or of 



Fig. 11. Same as Fig. QT5| but now the ordinate is 
log[K; nt / Qb[aLcr(L)] 2 ]. The solid lines ( Kmt x d 2 ' 11 ) have rj — 
and 77(n = -2) « 0.3, T]{n = 0) * 0.5, rj(n = +1) * 0.6 The 
case n = —2 shows also the same data furnished with scatter 
bars and shifted by a factor 10 3 for clarity. The scatter bars of 
the other two cases are of the same size. 



Kint, as argued in Sec.|5]and exemplified by our results for Kmt- 
Relevant for the dynamical evolution of these physical quan- 
tities is not the bulk velocity field, but the relative velocity 
(that is, the velocity gradient), which does not suffer this R— 
dependence. 
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Fig. 12. For n = +1, r n \ = R/16, cr(L) > 3: the uppermost data 
are scaled like in Fig. ^| the lowermost data like in Fig. 
For the other values of n, the behavior of the data is the same. 




-2-101234 



Fig. 13. Abscissa: log(l + 6). Ordinate: logA-int (arbitrary nor- 
malization) for the flat CDM model in a box of sidelength 
128 Mpc. For clarity, only some coarsening lengths are plot- 
ted (decreasing top down). The solid lines (/<i nt oc 5 2 ~ n ) have 
T] = and 77 = 0.5. 



In the cases n = 0, +1, the theoretical predictions and the 
scaling are well followed except at the earliest times and 
smallest cell sizes, when resolution effects are expected to be 
most important. An interesting result is the linear dependence 
of K mac with 6, Eq. Jl 3i . with a proportionality factor which 
according to Fig.[7]does not seem t o depe nd sensitively on the 
spectral index ra. lSeto & Sugivamal 1 1200 11) have studied K miC in 
the cases n — 0, +1 in the quasilinear regime (0.2 % cr % 1.2), 
which we have not addressed at all. 

The internal kinetic energy, K[ nt , is more sensitive to the 
small scale dynamics than K mac is; correspondingly, resolution 
effects are found to be more important than for K mac . However, 
only at the earliest time and smallest coarsening cell sizes for 
n — 0, +1 do they render the results unreliable. That's why the 
theoretical prediction for K mt {6) in the linear regime could not 
be tested when n = 0, +1; when n = — 2 the reason is that the 



Fig. 14. Abscissa: log(l + 6). Ordinate: log K mt (arbitrary nor- 
malization) for the flat CDM model in boxes of sidelength 
128 Mpc and 512 Mpc, respectively. In the later case, Km was 
multiplied by a factor 10 3 . For clarity, only some coarsening 
lengths are plotted. The solid lines (K; nt oc S 2 ^) have 77 = and 
77 = 0.5. 

linear regime was not really probed, which can be traced back 
to a finite-size effect. 

In the nonlinear regime, the function Ki nt (6) exhibits an 
interesting behavior. The virial prediction, Eq. l|8}, was ob- 
served only for n = 0, +1 asymptotically in the large-<S end 
of the curves. Otherwise, a polytropic-like dependence dl5> 
was found to fit better the data, with an exponent which does 
not seem to depend on time or coarsening length, only on the 
spectral index n. The same polytropic-like dependence is found 
for CDM models, albeit with a scale-dependent exponent , con- 
firming the results of an earlier work (Dommguez 2003). The 
values of the exponent 77 f or the CDM models analyzed here 
and by DommauezJ J2003I) are consistent with each other and 
with those of the self-similar models, in spite of the differ- 
ences in the number of parti cles in the simulations (and, com- 
pared to Dommguez 2003, the simulation algorithm itself). 
This suggests that the polytropic-like relation is not an arti- 
fact of the simulations, which in this respect seem to reproduce 
acceptably the limit N — > +00. Another hint in this direction 
is the simple relation which connects the results of the CDM 
models in boxes of different size, Sec. 14.41 in the larger box, 
R = 512 Mpc, a coarsening cell of a given mass contains less 
particles that in the smaller box, R = 128 Mpc. Nevertheless, 
this mass-resolution effect does not alter the functional de- 
pendence K[ m (6) at all, and can be accounted for by a scale- 
independent const ant offset. 

In the work by Kep ner et"al] Jl997l) . /Ci nt is also measured in 
CDM simulations (in their notation, K ml = gcrl/2). They find 
a polytropic-like dependen ce too , but w ith a slightly smaller 
exponent, 77 * 0.1 JKepner et ail 1 19971 Eqs. 17-18). We be- 
lieve this discrepancy to be a consequence of their simulation 
having too few particles (N = 32 3 ): as a consequence, they 
measured the function Ki nt (5) with coarsening cells having at 
most 100 particles (in one case); i n many cases, the c ells have 
less than a few tens of particles JKepner et alJI 19971 Fig. 3). 
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For comparison, the polytropic-like dependence in Fig. 1 1 31 is 
detected in coarsening cells containing a number of particles 
spanning ranges as wide as 30 - 7000 or 500 - 30000 (see also 
lDomfnguezl200 3). The work by Kee ner et alJJl997l) was moti- 
vated by comparison with redshift surveys. It relied on the cos- 
mic virial theorem and particular emphasis was put on the de- 
pendence with cosmological parameters (Cl m , Q.a). Our results 
show that departures from the virial predic tion are no t small at 
all, so that the method devised bv lKepner et al.N 19971) must be 
adjusted. More generally, our results warn against a straightfor- 
ward use of the cosmic virial theorem to estimate cosmological 
parameters from observations without first assessing that the 
employed observational data do indeed pertain virialized struc- 
tures. 

In the work by Nagam ine et alJ d200lb . the cosmic Mach 
number ( = K mac / Kmt in our notation, and K; nt = go 2 /2 in theirs) 
is measured in a ACDM hydrodynamical simulation, for three 
different length scales and as a function of the density. As a 
side-result, they also find a polytropic-like dependence for the 
velocity dispersion of groups of DM halos and galaxies (with 
77 « 0.5 - 0.7) — the authors do not elaborate much on this 
result. One must keep in mind that, compared to our simu- 
lations, theirs involves also the baryonic component and the 
formation of galaxies, which can affect the velocity dispersion 
jTissera & Domfnguez-Tenreiroll99"8T) . 

One can conceive two natural extreme cases of a 
polytropic-like dependence: the "virial" case, Tj^i^ - 0, when 
velocity dispersion is fixed by the local mass density, and the 
"isothermal" case, //isothermal = 1, when velocity dispersion is 
fixed by an external cause, e.g. tidal forces, free flow,... The 
values of 77 that we measure invariably fall between and 
1; the corresponding relations K mX (S) can be arguably under- 
stood as the outcome of the competition of the two effects 
("local virialization vs. global thermalization"), whose rela- 
tive strength varies with the spectral index n and the cell size 
and mass. However, an elaborated theory is required to lend 
support to this explanation, the ultimate goal being the "post- 
diction" of the relation Jl 5I >. T he value 77 = 1/3 was de- 
rived theoretically by Buc hert & Domfnguezl |l998), but we 
think this is irrelevant to our results, since certain restrictive 
assumptions were made (vanishingly small and isotropic ve- 
locity dispersion, approximately shear-free velocity field u), 
which are unlikely to hold in the regime where we find the 
polytropic-like dependence. The results from the simulations 
cannot be explained by any theory whose starting point is the 
usual thermodynamical theory or, more general ly, the (grand- 
canonical ensemble of statistical mechanics dSaslaw & Fang 
ll996HHochberg & Perez-Mercaderll996llde Vega et alJl9 98). 
since in that framework the kinetic energy is an extensive vari- 
able: /Cj nt = Tg (T is the kinetic tem perature) an d 77 = 1. As 
a side-remark, we notice that Sas law et al.l ( 1990) compute the 
velocity distribution allegedly in the framework of thermody- 
namics: but they use contradictory arguments and obtain in- 
stead that the kinetic energy scales like g 2 (that is, 77 = 0), and 
a velocity distribution different from the Maxwellian one char- 
acteristic of thermal equilibrium and which should follow from 
the (grand-)canonical ensemble probability. 



The discovered relationship K mt (g) is useful for an im- 
prov ed model of structure forma tion by gravitational instabil- 
ity dBuchert & Domfnguezll99'8T) : the dust model (pressureless 
fluid) is added a term proportional to the gradient of K; nt (a ki- 
netic pressure), in order to account for the reaction of the dy- 
namically generated velocity dispersion on the evolution. The 
evolution equation for the velocity field u(x, t) then reads 

<9u a 1 1 

_ = — u - -(u V)u + w V#£i nt , (16) 

at a a ag 

where the peculiar gravitational acceleration w is given by 
Poiss on's equation. Fur t her theoretical studies of this m odel 
dAdler & Buchertl 1 19991: IBuchert et alJ 1 19991: IMaartens et alJ 
Il999t iMorita & Tatekawall200ll ITatekawa et alJl2002l) work 
with a K mt which depends only on density, e.g. a pure 
polytropic-like dependence with values of the exponent 77 in 
concordance with our measurements. Our results show that the 
functional dependence of K m on g is somewhat more com- 
plicated than purely polytropic and changes with time and 
coarsening length. Nevertheless, since the term -V/f; nt is a 
pressure, it opposes compression in collapsing regions. More 
can be learned about the beha vior of this term when some 
simplifications are introduced (IBuchert & Domfnguezl Il998l 
IBuchert et alJll999l) : one assumes that the evolution follows 
the dust model prediction in the form of the Zel'dovich ap- 
proximation (basically that u oc w) "almost everywhere", i.e., 
except near potential density singularities where the effect of 
— V/fint becomes relevant. One can then apply boundary layer 
theory to show dDommguezll200ol) that this term does behave 
"adhesively" and indeed succeeds in preventing the formation 
of a singularity, provided Ar; nt is a function of g and K mt (g)lg is a 
growing function of g (meaning 77 < 1 for a polytropic-like de- 
pendence). This behavior is robust against the observed time- 
dependence in the relation K mt (g), being much slower than the 
time-scale of collapse. One can conclude that the dependence 
Kintfe measured in N-body simulations leads to the same 
qualitative "adhesive" behavior as the simpler dependences ad- 
dressed theoretically in the literature. 

When it comes to inserting our results in the theoretical 
model d!6i . there are some issues which we have not addressed 
but may be relevant to a better understanding of the model. 
First, it must be noticed that the average relationship K ml (g) 
does not mean in principle a one-to-one dependence between 
K ml and g; on the contrary, the data scatter around the average 
dependence, Fig.^] In fact, the derivation of Eq. ( II6J 1 yields 
in reality a term -VK mt (Buchert & Domfnguez 1998): the in- 
fluence of the scatter on the model outputs should be quanti- 
fied and, if pro ven relevant, incorporated in the model, e.g. as 
a noisy source (Buchert et alJll999l) . Another issue of possible 
concern is the amount of velocity dispersion in the coarsen- 
ing cells associated to "bound structures" (as opposed to the 
amount associated to particle flow between neighboring cells); 
in this context, it would also be interesting to assess the contri- 
bution to velocity dispersion from "ordered motion", e.g. due 
to a net angular momentum. 

In conclusion, we have studied the density dependence of 
the macroscopic and internal kinetic energies in coarsening 
cells. We could identify the influence of finite-size and reso- 
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lution effects on the measured physical quantities. When these 
effects were irrelevant, we could confirm some of the theoret- 
ical asymptotic predictions. Finally, we found that in an inter- 
mediate range of densities, the velocity dispersion scales as a 
power of the mass density, with an exponent different from the 
virial prediction. 
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Appendix A: Estimates with the linear solution 



x[W(Lk)W(Lq) - W(L(k + q))] <J mic (k, t)6 nic (q, t). 



In this Appendix we collect the mathematical calculations 
which lead to Eqs. J4I6I . The main idea is that K lnt and /c mac 
are determined by the dominant contribution of modes in the 
linear regime when L » r„i, so that they can be estimated by 
inserting the linear solution in the definitions Q. These defini- 
tions can be rewritten as follows: 

g(x,t) = J" ^ W {^-)g mic (y,t), 

/dy / x — y \ 
U W [~ ~ \Qmic^mic{y,t), 

K mt {x,t) w(^)e«fc«L(y.O- ^gu 2 (x,t), 

in terms of the formal microscopic fields 4 
g mic (x, t) := —3 V S°\x - x a (t)), 

m v-i m 
g mic ii mic (x,t) := —3 2^^(05 (x-x„(f)). 

v / a 

We introduce the Fourier transform of any spatial field <f>(x), 
denoted by a tilde and defined as 

0(k):= Jdxe' kx <^(x). 

The velocity then reads 



gu(x, t) 



dq , 



e-'' (k+q) ' x g mic (k, f)u m/c (q, t)W(L(k + q)), 



and the internal kinetic energy 

*m(*,0 = J f ^?e-" h «)-6„ c (k,0x (A.1) 
2 J (27r) y 

x[u m/c (q, f) • u, rac (p, t)W(L(k + q + p)) - u(q, t) ■ u(p, t)W(Lk)]. 

Even if L <c r„i, these expressions contain contributions from 
nonlinear modes (r n \k » 1). The hypothesis (Sec. [5J is that 
these contributions are nevertheless negligible compared to 
those from the linear modes (r a \k <sc 1). This implies u(k) « 
u„„ c (k)W(Lk) in the linear regime, and the growing linear so- 
lution in an Einstein-de Sitter background yields in turn 



ik- 

u m , f (k,f) ~ -a-j^6, nic (k, ?), 



(r nX k «: 1) 



with 6 m ic — (gmiclgb) - 1- Inserting these linear relations in 
the definition JA. II we obtain, to lowest order in the inhomo- 
geneities, 



#int(x, 



1 2 fdkdq 

2 ebh J W 



k 2 q 2 



X 



4 In the mathematical manipulations in this Appendix, one can as- 
sume that the Dirac delta is regularized by a short-distance cutoff 
which is set to zero in the final results. In this way, the field iw(x) 
is defined, and not only the product Q mjc u mic (x). 



For Gaussian initial conditions, the average of Eq. Q can 
be easily computed using standard techniques for the Gaussian 
functional integrals (e.g. Zinn- Justinll996l) : 



(6 (3 \g(x,t;L)-g)} 



■ exp 



yl2ng b {i)cr{L,i) \ 2cr 2 {L,t), 
(S mic (k, t)6 mic (q, t)6 (3) (g(x, t; L) - g)) 



{6<V(g(x,t;L)-Q)) 



(2nfP(k, ?)<5 (3) (k + q)- 



—27j-Z f 1 - ^7T^) P ^ OW(ik)W(Lq)e' (k+ i ) x . 
cr z (L, t) \ cr l (L, t) ] 

With these expressions, one recovers the result ijfji with the co- 
efficients given by 



<A' inl )(L.n = -g b a 2 J ^ [1 - \W(Lk)\\ (A.2) 



B{L, t) = 



gba 



f dkc 



dq k q 



2cr 2 (K int )(L, t) J (2nfk 2 q 2 
xP(k, t)P(q, t)W(Lk)W(Lq)W(L(k + q)). 



(A.3) 



These integrals are IR-convergent provided P(k — > 0) ~ k", 
-3 < «, since W(k — » 0) ~ 1 +o(k 2 ) (<-» the window function is 
normalized to unity and decays fast enough at large distances). 
On the other hand, there is an implicit ultraviolet (UV) cutoff, 
k,q < r nl ', because these expressions have been derived for 
the linear regime. The limit "r n i — > 0" of the coefficient B is 
finite, since we assume that W(k — > 00) — > fast enough (<-> 
the window is smooth enough: the decay is exponential for a 
Gaussian window, algebraic (~ Ar 2 cosA:) for a spherical top- 
hat window). (K- mt ) is however UV-convergent only if P(k — > 
00) ~ k", n < -1, indicating that it is determined by nonlinear 
modes when n > - 1 . 

A similar reasoning can be repeated for the macroscopic ki- 
netic energy: inserting the linear relationships in the definition 
-Kmac = (l/2)£>u 2 one gets 



1 , f dk 

I 6 '" J ^ 



dk m Wk)i 2 = <* mac >, 



) 3 k 2 



which is IR-divergent in the range of spectral indices -3 < n < 
-1. 



